home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 051-075 / disk_065 / prep / ktest.p < prev    next >
Text File  |  1992-05-06  |  4KB  |  172 lines
  1. c test program for Kutta.  Probably needs double precision.  If so
  2. c do:
  3. c    prep -d DOUBLE ktest
  4.  
  5. #include "premac.h"
  6. : DIM    2 ;
  7.  
  8.     implicit real ( a-h, o-z )
  9.     real y(DIM)
  10.     external simple
  11.  
  12.     xi = 10
  13.     xf = 0
  14.     y(1) = exp(xi)
  15.     y(2) = exp(xi)
  16.     n = DIM
  17.     del = .1
  18.     accurc = 1.e-8
  19.     imax = 10
  20.     
  21.     write(*,*) kutta( xi, xf, y, n, del, accurc, imax, simple ),
  22.      *        ' iterations'
  23.  
  24.     do i = 1, DIM
  25.        write(*,100) i, xf, y(i)-1
  26.     end do
  27.  
  28.     stop
  29. 100    format( ' error in y(', i2, ') at x=', g12.5, ' is ', g12.5 )
  30.     end
  31.  
  32.  
  33. c functions to be integrated
  34.     subroutine simple( x, y, f )
  35.     implicit real ( a-h, o-z )
  36.     real y(DIM), f(DIM)
  37.  
  38.     f(1) = y(2)
  39.     f(2) = 2*y(2) - y(1)
  40.  
  41.     return
  42.     end
  43.  
  44.  
  45.  
  46.  
  47. c Runge Kutta ODE solver, real version
  48. c
  49. c This routine solves a system of 1st order ODE's of the form
  50. c
  51. c     dYn
  52. c     --- = Fn(x,Yn)
  53. c     dx
  54. c
  55. c where the functions Fn are defined by the external routine
  56. c equa( x, Y, F ).  n must be <= DIM, defined below.
  57. c
  58. c variable defs
  59. c    xi, xf        integrate from xi to xf
  60. c    y        initial y's,  solutions also returned here
  61. c    n        number of initial y's
  62. c    del        initial dx
  63. c    accurc        accuracy measure
  64. c    imax        maximum times that del may be halved
  65. c    equa        external functions to be integrated
  66. c
  67. c returns the number of iterations, -1 on error
  68. c
  69. c originally written ??? by ???
  70. c converted to a readable form 2/2/87, P.R.Ove
  71.  
  72. #include "premac.h"
  73.  
  74. : INCREASING_DONE    ( (xn+h>=xf) & (xf>=xi) ) ;
  75. : DECREASING_DONE    ( (xn+h<=xf) & (xf<=xi) ) ;
  76. : DONE            ( INCREASING_DONE | DECREASING_DONE ) ;
  77. : DIM    6 ;
  78.  
  79. do limits [ n ]
  80.  
  81.     integer function kutta( xi, xf, y, n, del, accurc, imax, equa )
  82.     implicit real (a-h,o-z)
  83.     real    y(DIM), yi(DIM), yn(DIM),
  84.      *        k1(DIM), k2(DIM), k3(DIM), k4(DIM), k5(DIM),
  85.      *        f(DIM), e(DIM), f1(DIM)
  86.     real    test, xi, xf, del, accurc, xn, h
  87.     real    amax1, abs
  88.     logical quit
  89.  
  90. c initialization
  91.     if ( n > DIM ) then
  92.        write(*,*) 'KUTTA: too many equations: ', n, '>', DIM
  93.        kutta = -1
  94.        return
  95.     end if
  96.  
  97.     kutta = 0
  98.     n_halves = 0
  99.     xn = xi
  100.     h = del
  101.     quit = FALSE
  102.     yn(#) = y(#)
  103.  
  104. c fix sign of h if not sensible
  105.     if (  ( xf > xi & h < 0 )  |  ( xf < xi & h > 0 )  ) h = -h
  106.  
  107. c check to see if finished. adjust h so that y(xf) is returned.
  108.     begin
  109.        if ( DONE ) then
  110.           del = h
  111.           h = xf - xn
  112.           quit = TRUE
  113.        end if
  114.  
  115. c runge kutta calculation. calculates y(xn + h) from y(xn) and dy(xn).
  116.        call equa(xn, yn, f1)
  117.  
  118.        begin
  119.  
  120. [          k1(#) = h*f1(#)/3.
  121.           yi(#) = yn(#) + k1(#)    ]
  122.  
  123.           call equa(xn + h/3., yi, f)
  124.  
  125. [          k2(#) = h*f(#)/3.
  126.           yi(#) = yn(#) + k1(#)/2. + k2(#)/2.    ]
  127.  
  128.           call equa(xn + h/3., yi, f)
  129.  
  130. [          k3(#) = h*f(#)/3.
  131.           yi(#) = yn(#) + 3.*k1(#)/8. + 9.*k3(#)/8.    ]
  132.  
  133.           call equa(xn + h/2., yi, f)
  134.  
  135. [          k4(#) = h*f(#)/3.
  136.           yi(#) = yn(#) + 3.*k1(#)/2. - 9.*k3(#)/2. + 6.*k4(#)    ]
  137.  
  138.           call equa(xn + h, yi, f)
  139.  
  140. c accuracy check. repeats with h halved if check fails.
  141.           test = 0.0
  142.  
  143. [          k5(#) = h*f(#)/3.
  144.           e(#) = (k1(#) - 9.*k3(#)/2. + 4.*k4(#) - k5(#)/2.)/5.
  145.           test = amax1(test, abs(e(#)))    ]
  146.  
  147.        while ( test >= accurc )
  148.           n_halves = n_halves + 1
  149.           if( n_halves >= imax ) then
  150.              del = h
  151.              kutta = -1
  152.              write(*,*) 'KUTTA: step made too small:'
  153.              write(*,*) '       del = ', del, '    at x = ', xn
  154.              return
  155.           end if
  156.           h = h/2.
  157.           quit = FALSE
  158.        again
  159.  
  160. c increase h to 2h if acceptable. return if finished.
  161.        yn(#) = yn(#) + (k1(#) + 4.*k4(#) + k5(#))/2.
  162.        xn = xn + h
  163.        kutta = kutta + 1
  164.        if ( test < accurc/32. ) h = 2.*h
  165.  
  166.     until ( quit )
  167.  
  168.     y(#) = yn(#)
  169.  
  170.     return
  171.     end
  172.